###############################################
#FIGURE 1
###############################################

library(estimatr)
library(ggplot2)
library(grid)

####load cleaned data####
data1 <- read.csv("./data/data1.csv")
data2 <- read.csv("./data/data2.csv")
data3 <- read.csv("./data/data3.csv")


####results without controls####

lm1 <- lm_robust(prej_direct_1 ~ as.factor(trt), data=data1)
lm2 <- lm_robust(prej_direct_1 ~ as.factor(trt), data=data2)
lm3 <- lm_robust(prej_direct_1 ~ as.factor(trt), data=data3)


com1 <- c("Identity Threat\nvs. Control", lm1$coefficients[2], (lm1$coefficients[2] + (1.65*lm1$std.error[2])), (lm1$coefficients[2] - (1.65*lm1$std.error[2])),
          (lm1$coefficients[2] + (1.96*lm1$std.error[2])), (lm1$coefficients[2] - (1.96*lm1$std.error[2])))
com2 <- c("Whataboutism\nvs. Identity Threat", lm2$coefficients[2], (lm2$coefficients[2] + (1.65*lm2$std.error[2])), (lm2$coefficients[2] - (1.65*lm2$std.error[2])),
          (lm2$coefficients[2] + (1.96*lm2$std.error[2])), (lm2$coefficients[2] - (1.96*lm2$std.error[2]))) 
com3 <- c("Whataboutism\nvs. Control", lm3$coefficients[2], (lm3$coefficients[2] + (1.65*lm3$std.error[2])), (lm3$coefficients[2] - (1.65*lm3$std.error[2])),
          (lm3$coefficients[2] + (1.96*lm3$std.error[2])), (lm3$coefficients[2] - (1.96*lm3$std.error[2]))) 
com <- data.frame(rbind(com1, com2, com3))
colnames(com) <- c("treatment", "effect", "cih90", "cil90", "cih95", "cil95")

com$effect <- as.numeric(com$effect)
com$cih90 <- as.numeric(com$cih90)
com$cil90 <- as.numeric(com$cil90)
com$cih95 <- as.numeric(com$cih95)
com$cil95 <- as.numeric(com$cil95)
com1a <- com


####results with controls####

lm1 <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1)
lm2 <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm3 <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3)


com1 <- c("Identity Threat\nvs. Control", lm1$coefficients[2], (lm1$coefficients[2] + (1.65*lm1$std.error[2])), (lm1$coefficients[2] - (1.65*lm1$std.error[2])),
          (lm1$coefficients[2] + (1.96*lm1$std.error[2])), (lm1$coefficients[2] - (1.96*lm1$std.error[2])))
com2 <- c("Whataboutism\nvs. Identity Threat", lm2$coefficients[2], (lm2$coefficients[2] + (1.65*lm2$std.error[2])), (lm2$coefficients[2] - (1.65*lm2$std.error[2])),
          (lm2$coefficients[2] + (1.96*lm2$std.error[2])), (lm2$coefficients[2] - (1.96*lm2$std.error[2]))) 
com3 <- c("Whataboutism\nvs. Control", lm3$coefficients[2], (lm3$coefficients[2] + (1.65*lm3$std.error[2])), (lm3$coefficients[2] - (1.65*lm3$std.error[2])),
          (lm3$coefficients[2] + (1.96*lm3$std.error[2])), (lm3$coefficients[2] - (1.96*lm3$std.error[2]))) 
com <- data.frame(rbind(com1, com2, com3))
colnames(com) <- c("treatment", "effect", "cih90", "cil90", "cih95", "cil95")

com$effect <- as.numeric(com$effect)
com$cih90 <- as.numeric(com$cih90)
com$cil90 <- as.numeric(com$cil90)
com$cih95 <- as.numeric(com$cih95)
com$cil95 <- as.numeric(com$cil95)
com1b <- com

####combined figure####

com1a$Controls <- "Without"
com1b$Controls <- "With"
full1 <- rbind(com1a, com1b)

level_order <- c("Identity Threat\nvs. Control", "Whataboutism\nvs. Identity Threat", "Whataboutism\nvs. Control") 

ggplot(full1, aes(x = factor(treatment, level = level_order), y = effect, color=Controls)) +
  geom_point(stat="identity", size=4, position = position_dodge(width = 0.15)) + 
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width = 0.15,),
                 linewidth = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width = 0.15),
                 linewidth = 2) +
  ylab("Treatment Effect") +
  xlab("Treatment Type") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  theme_bw()+
  theme(
    axis.title=element_text(size=18), 
    axis.title.x = element_text(margin = unit(c(7, 0, 0, 0), "mm")),
    axis.title.y = element_text(margin = unit(c(0, 7, 0, 0), "mm")),
    plot.title = element_text(size=26),
    axis.text=element_text(size=16)
  )    

ggsave("fig1.png", width = 10, height = 6, dpi = 300)

####additional statistics related to these results reported in-text####

#control group mean
data_clean <- read.csv("./data/data_clean.csv")
datac <- subset(data_clean, data$trt == "control")
mean(datac$prej_direct_1, na.rm=T)

#whataboutism vs identity threat effect size in SDs
lm2$coefficients[2]/sd(datac$prej_direct_1)

#whataboutism vs control effect for over-40s
data3_over40 <- subset(data3, data3$age>40)
summary(lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_over40))

#whataboutism vs control effect for above-average pre-treatment national id
data3_highnid <- subset(data3, data3$natid_covar_pca>mean(data3$natid_covar_pca, na.rm=T))
summary(lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_highnid))
